Gaussian process forecasts Pseudogymnoascus destructans will cover coterminous United States by 2030

Abstract White‐nose syndrome has been decimating populations of several bat species since its first occurrence in the Northeastern United States in the winter 2006–2007. The spread of the disease has been monitored across the continent through the collaboration of many organizations. Inferring the rate of spread of the disease and predicting its arrival at new locations is critical when assessing the current and predicting the future status and trends of bat species. We developed a model of disease spread that simultaneously achieves high‐predictive performance, computational efficiency, and interpretability. We modeled white‐nose syndrome spread using Gaussian process variations to infer the spread rate of the disease front, identify areas of anomalous time of arrival, and provide future forecasts of the expected time of arrival throughout North America. Cross‐validation of model predictive performance identified a stationary Gaussian process without an additional residual error process as the best‐supported model. Results indicated that white‐nose syndrome is likely to spread throughout the entire continental United States by 2030. These annually updatable model predictions will be useful in determining the horizon over which disease management actions must take place as well as in status and trend assessments of disease‐affected bats.

management and intervention strategies. However, spatial models can be complex and computationally intensive, making it difficult to reproduce, to update over time, and to share among practitioners.

White-nose syndrome (WNS) is an emerging infectious disease
in North American hibernating bats that has driven widespread population declines and currently threatens several species with local or global extinction (Cheng et al., 2021;Frick et al., 2010;McGuire et al., 2016;Thogmartin et al., 2013;Thogmartin, King, McKann, et al., 2012;White-nose Syndrome Response Team, 2021a).  et al., 2015) and visualized at https://www.white noses yndro me.org (White-nose Syndrome Response Team, 2020), demonstrates the state of the disease front and its diffusive spread outward from its initial detection (Figure 1).
When WNS arrives at a bat hibernaculum, in most cases the winter survival rate of the population plummets for several consecutive years (Cheng et al., 2021). Based on observed spread over the past decade and a half, it seems likely that the fungal pathogen will continue to spread over the entirety of the United States and possibly the entire continent. Uncertainty remains, however, because although many known sites have been sampled for WNS, many more remain unsampled, in addition to sites still remaining to be discovered. When assessing the status and trend of WNS-susceptible species, knowing where these populations are in the disease invasion process is crucial . For those sites for which abundance monitoring data exist but the presence of WNS or Pd has not been determined, these predictions are essential for correctly describing population status. These data can inform models as a covariate in changepoint analyses or for parameterizing vital rates within demographic models (Erickson et al., 2014;. NABat aims to provide a periodic assessment of species status and trend, and for cave-hibernating species subject to the deleterious effects of this disease, predictions of disease arrival, with uncertainty, will be crucial for credible assessments. Thus, the primary goal of this work is to produce spatiotemporal predictions of the year of disease arrival within North America, predictions that are readily updatable as new surveillance data become available. Several studies have attempted to model the spread of Pd/ WNS across space and time. Various approaches have used statistical models to understand the fundamental processes driving the disease and its effect on animal populations. Inference under these F I G U R E 1 The observation locations, determination status, and the observed year of arrival of WNS or Pd. the Pd determination status indicates presence of the pathogen whereas WNS determinations were made with a histological examination. Raw data available at https:// www.white noses yndro me.org. models has yielded estimates of the origin, growth, and spread of the disease, ecological characteristics predicting spatiotemporal disease arrival, and the mortality and behavioral effects it has caused in bat populations. On the continental scale,  use a spatiotemporal, mixed-effects logistic regression to predict WNS arrival, whereas Maher et al. (2012) andKramer et al. (2021) modeled the spatial probability of infection using logistic regression gravity models. Meierhofer et al. (2021) investigated the spread of WNS in Texas using host-pathogen susceptible-exposedinfected (SEI) differential equation models. Wilder et al. (2015) modeled risk of infection among little brown bat Myotis lucifugus hibernacula using both geographic proximity and genetic similarity, finding genetic distance improved model fit. Also using genetic data, Thapa et al. (2021) Hefley et al. (2017Hefley et al. ( , 2020 modeled the growth and diffusion of the pathogen using partial differential equations, which provide a convenient method for incorporating science-based physical models, but numerical algorithms discretizing these equations and performing estimation of parameters can be computationally expensive for even simple spatial models (Hooten et al., 2013;Hooten & Wikle, 2008;Wikle, 2003).
Physical diffusion processes are inherently stochastic. Ecological diffusion models incorporate this physical model, but often include additional elements of uncertainty. For example, the bats most affected by WNS are migratory, but the details about seasonal migrations, the mixing of individuals, and the interaction among species and resulting interspecific spread are not well understood; neither are environmental factors that might mediate spread such as the type of hibernaculum or roosting site (from traditional caves and trees to human-made structures such as culverts) or the medium of transmission, including the possibility of human-induced spread.
Thus, phenomenological models are useful, particularly when one goal of the analysis is predictive performance. Grider et al. (2021) provides an example prioritizing prediction, using a generalized additive model in the form of a logistic regression on presence of Pd including separable, smooth functions of space and time. Results indicate WNS spread is best predicted when including smoothed terms of space and time.
These previous modeling efforts have provided useful information adding to our knowledge of WNS but are not adequate for conservation planning because they are not easily updateable or repeatable. Repeatability can be limited by data availability (e.g. models incorporating logistic regression require presence/absence data and SEI models require host-pathogen data) or by model complexity (e.g., growth and diffusion models require extensive computational resources and statistical expertise). In particular, Maher et al. (2012) and Meierhofer et al. (2021) mention the importance of cave density and environment. However, this information is not publicly available.
In addition, the locations of hibernacula in the western United States are not as well known or documented, limiting the predictive power of such models in the area where predictions are most needed. The uncorrelated Gaussian error structure of these models rely only on covariates to capture the variation in the spread and the anomalous observations such as in Washington, which is often inadequate for predictive purposes.
In the absence of information about cave density and environment or genetic data, a distance-based model can account for anomalous detections using a generalized error structure. A Gaussian process (GP) can achieve this by including a correlated error structure, modeling local anomalies and subsequent spread, and also possibly by including a residual error structure. The GP is flexible and accurately captures non-diffusive spread, whereas human-mediated transmission (e.g., Washington) would not be explained using more covariates such as cave density or environment.
In this study, we use publicly available data reporting year of Pd arrival or WNS emergence to develop several GP models that are fast, reproducible, and easily updated with each new iteration of data (U.S. Geological Survey, 2021; White-nose Syndrome Response Team, 2020). From these models, we obtain a predicted year of Pd/ WNS arrival over space, as well as inference about the rate of disease spread. Given the rapid spread of Pd/WNS and associated effects to bat species, it is critical to provide predictions that can be updated annually at the same pace of the WNS disease cycle and that can be easily reproduced and shared among conservation practitioners.
Although our model is specific to the WNS, the principles underlying our modeling effort are translatable to other disease systems facing similar conservation pressures. Indeed, GPs have been used in epidemiology to capture uncertainties in both the measured data and the underlying disease spread process (Ak et al., 2018;Bhatt et al., 2017;Vanhatalo & Vehtari, 2007). GPs can include ecological and environmental covariates while accounting for additional spatially autocorrelated variation, which can result from either omitted covariates or because ecological processes are often inherently spatial (Wright et al., 2021). Thus, GP models are well-suited to the Pd/ WNS monitoring data that are publicly available.
We first present the general model and the variations we implement in Sect 2. We then compare inference and prediction results of the model variations in Sect 3, present a retrospective crossvalidation experiment assessing the predictive accuracy of each model as observations are made over time, and provide details on computational performance. We conclude with a discussion of the significance, limitations, and possible extensions of these results in The four statistical models we compare in this work can be written as variations of the following process model: where (s) ∼  0, K s, s � and z(s) ∼  0, K z s, s � are mean zero GPs with covariance functions K ε and K z . A review of GP estimation and prediction can be found in Gelfand et al. (2010), Cressie (1993), and Rasmussen and Williams (2005).
The first of four models that serves as a baseline comparison is the linear model, which is obtained by setting z = 0 and white noise covariance K s, s � = 2 ⋅ s = s � , where is the indicator function.
Using vector notation we can write the data distribution of obser- parameters to be estimated are β and the residual variance τ 2 .
The mean of the linear model is [y(s)] = + x(s) . μ is included as an offset parameter, fixed at the year 2007 for all observations, which implies that β represents an unknown spread rate parameter to be estimated in units of years per kilometer (km). In the results section, we report the inverse 1∕β in the more familiar units of km/year. This combination of fixed effects was chosen to preserve the interpretation of the effect of distance from Howes Cave as a constant radial spread rate since the introduction of the disease in North America. The model mean is the same for the three GP models about to be described; complexity is added through the secondorder structure of the processes.
The second model, the stationary GP, does not include ε but instead includes z, a spatially autocorrelated, stationary, mean zero GP with the Matérn covariance function where d is the great circle distance between two spatial locations of interest s and s ′ . The data flexibility among determination statuses via different variance parameters may indicate whether samples taken in proximity but with different determination statuses (or sampling methods) are likely to agree.
In all models, temporal dynamics are modeled in the process mean by the linear predictor defined with μ and the radial spread rate β, and the remaining variance is partitioned into mean zero processes ε and z endowed with specific second order structure.
Deviations from this radial spread are captured by the spatially correlated GP z and the residual error process ε. Including a residual error process implies the belief that there is either uncaptured mi- Likelihood methods yield estimates of the fixed effect β and covariance parameters, and predictions with uncertainty bounds can be obtained at arbitrary spatial locations for all models. All parameters were estimated except the smoothness parameter in the three GP models, which was fixed at ν = 0.5. In R (4.1.0), we use lm to fit the linear model and the optim function to minimize the negative log restricted likelihood function of the stationary, nugget, and heteroskedastic GPs (R Core Team, 2021).
The primary objective of this analysis is to provide the most ac- were obtained using the best linear unbiased predictor for spatial data (Cressie, 1993), the method otherwise known as regressionkriging or kriging with external drift.

| Status and forecast of WNS
The GP models indicate the disease could cover the coterminous United States by 2030.

The parameter estimates obtained for the linear model and three
GP models using all data available are given in Table 1 The four models we examined provided increasingly complex insight into disease and fungal spread ( Table 1). The linear model indicated the rate of spread was fairly precise at about 182 km/year. The estimate from the GP model without nugget effect was nearly the same but less precise. Conversely, the nugget and heteroskedastic GP models indicated a much more rapid rate of spread (207 km/yr), but prediction intervals were much wider.
The models differed in their estimation of parameters defining the residual error process ε. The linear model suggested disease and fungal arrival may have been more than two years earlier or later than observed based on the estimate of τ. The residual error estimate for the nugget GP model was smaller, which might be expected as some variance in arrival year in the monitoring observations can be partitioned to the autocorrelated process standard deviation σ.
The nugget GP estimate of τ indicated slightly more than one year residual error, while the heteroskedastic GP model estimates of τ 1 ,τ 2 , and τ 3 indicated arrival may have occurred 2.41 years earlier or later than observed at Pd-positive sites, 0.61 year at WNS-positive sites, and 1.19 years at WNS-suspect sites (see Table 1).
Estimates for the spatial autocorrelation parameters σ and θ have a clear dichotomy between the stationary GP model and the nugget and heteroskedastic GP models. Estimates for σ are 15%-30% larger for the nugget and heteroskedastic GPs compared to the stationary GP. The estimated correlation ranges showed an even greater contrast: the estimate for the GP without nugget was 79 km whereas the estimates for the nugget and heteroskedastic GP models exceeded 900 km.
For the stationary GP model, we calculated mean predictions and standard errors on a finely spaced grid for visualization ( Figure 2).
The latent process, excluding the residual error variation, represents the estimated year of arrival of Pd/WNS at a given spatial location.
See the following section for justification of the choice to display the stationary GP model results. Also see in the supplementary material Figure S3 for predictions from all four models as well as Figure S4 for model-weighted predictions. Note that plots of surface predic-

| Cross-validation
To assess the retrospective cross-validation of our models, for each year from 2012-2021 we train the models on all data up to the given year and then perform pseudo out-of-sample forecasting on all the data collected the year following the given year.
Estimates of the spread rate for the four models converged on 150-350 km/year (with uncertainty bounds) for all four models by 2021 (top row of panels in Figure 3). The middle row of Figure 3 shows standard deviation parameter estimates over time (in units of years), with τ in the first column, σ in the second column, and τ 1 , τ 2 , and τ 3 in the third column. The bottom row of Figure  This cross-validation experiment reveals how model estimates evolve over time as new data are collected. Model weights calculated from the CRPS predictive performance metric indicated which models might be useful for making predictions. Model averaging using these weights provides one way to combine models with several hypotheses into probabilistic predictions (Nichols et al., 2021). Figure S3 shows predictions using all available data from 2007-2022 for the four models we investigated, and Figure S4 combines these predictions weighted by the 2021 cross-validation weights.
All stationary GP models took under a minute to fit, and all heteroskedastic GP models took less than 2.5 minutes to fit on a Dell Precision 7550 with 10th Generation Intel Core i9-10885H vPro processor. Predictions on a fine spatial grid with over 180,000 spatial locations took about 1.5 minutes for all models, and prediction intervals on the same grid took under 45 minutes using 16 cores for all GP models fitted using monitoring data from all years.

| DISCUSS ION
In this study, we investigated several Gaussian process models pro- America. Programmatically, this work provides capacity for NABat (Loeb et al., 2015) with predictions incorporated into annual status and trend assessments of bat species; as new observations become available, the speed and ease of calculation allow rapid incorporation into assessments. Quantitative results from models are a key component of structured decision-making and adaptive resource management (Nichols et al., 2021;Runge et al., 2020). Being able to predict disease arrival could also help inform preventative measures to protect bat populations (Grider et al., 2021;Meierhofer et al., 2021). These results benefit policymakers, researchers, and the public by providing a better understanding of what the future outlook for bats may look like, and what conservation actions may need to be taken to protect these species (Bernard et al., 2020).
These approaches, however, can be hampered by a paucity of infor- The second-order model components of the GP models can also be interpreted within the epidemiological context. The residual error, for example, could represent measurement error or microscale variation from missing covariates, and regions with spatially correlated deviations from the estimated spread rate are attributable to anomalous early or late disease arrival. Another advantage of this simple approach is that the only pieces of information needed to predict to unobserved locations are the spatial coordinates and the distance from the disease epicenter, which is easily calculated and time invariant. Our exclusion of additional covariates at this time is intentional for two reasons: predictions requiring these covariates may limit their spatial or temporal support, and the meaning of spread rate within the model can become obscured, especially with the inclusion of nonlinear functions of space and time.
However, the autocorrelation process z in the GP models is, however, able to account for deviations from a constant spread rate without additional covariates. The slow spread estimates of the GP model seems to be accurately capturing the observed spread exhibited by the data. This is a boon for using a GP model to include a spatial autocorrelation process in comparison to competing regression-based models based on covariates alone without a generalized error structures. The GP models are able to incorporate multiple components of variation existing within the monitoring data through the autocorrelated and uncorrelated processes z and ε. The different estimates of the extent of spatial autocorrelation in the data measured by σ and θ are illuminated by comparing the stationary GP with the nugget and heteroskedastic GPs. In the heteroskedastic model, differences in residual standard deviation estimates indicate how data from each of the three determination statuses contributed to the spatial identification of the disease front.
These results tell a compelling story regarding the spread of Pd/WNS in North America, identifying areas significantly deviating from the concentric spread such as the late arrival in the Southeast and early arrival in Washington, where the spread may have been human-mediated. As the spread of the disease becomes more erratic over time, however, anomalies such as that on the eastern coastal plain may not be completely captured by the model. One possible explanation is that testing improved over time, allowing Pd to be detected sooner, but during the initial years of model fitting, The choice of model used for disease spread ultimately depends on the required output. The GP models we present can be readily fitted on an annual basis as new data become available. However, if fine-grained spatial information is required about fungal growth or diffusion, an ecological diffusion model may be preferable.
The spatiotemporal prediction forecasts we produced are valuable to those studying bats in areas where Pd/WNS has yet to arrive but may have a significant effect in the future, such as the southwestern United States. Several features that distinguish the ecological niche of bats in the western United States make it difficult to make forecasts in this region: there may be fewer WNS-susceptible species and bats hibernate at lower densities often in unconventional, inaccessible, or unknown sites (Bogan et al., 2003;Hendricks, 2012). In addition, the dynamics of the spread from the long-distance transmission to Washington are still unfolding and contributes to uncertainty in the disease front in the West. In efforts not shown, we explored one model variation including a distance covariate from the first observation in Washington, but this model variation did not appreciably improve model fit or performance and thus was not explored further; in future years, however, this model remains a viable candidate for further examination. writing -review and editing (equal).

A PPE N D I X A
Model estimates for standard deviation parameters shown in the bottom row of Figure 3 generally increase over time as more observations over a larger spatial domain are used to fit the models. To investigate the effect of the increasing area of the spatial domain on the magnitude of the standard deviation parameters, we plot the ratio of each standard deviation parameter to the area of the spatial convex hull of observations in Figure S1. These area weighted standard deviation estimates generally decrease initially before stabilizing over the final 5 to 10 years. Stable estimates may indicate these errors arise from irreducible background variability emerging from the monitoring data. If this is the case, more observations will increase the predictive performance of this approach but will not diminish the size of the irreducible errors. Stable estimates of standard deviation parameters over time similarly yields stable uncertainty quantification about model predictions, based on the monitoring tools, datasets, and the methods we used to conceptualize and model the problem. Figure S2 shows the number of determination status observations in the dataset each year. Note the first appearance of Pseudogymnoascus destructans (Pd) determination data in 2011 and the marked increase in 2017 and the following years.